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ABSTRACT 


The primary objective of this work was to demonstrate the feasibility of a new potential/viscous 
flow coupling procedure for reducing computational effort while maintaining solution accuracy. This 
closed-loop, overlapped, velocity-coupling concept has been developed in a new two-dimensional code, 
ZAP2D (Zonal Aerodynamics Program-2D), a three-dimensional code for wing analysis, ZAP3D (Zonal 
Aerodynamics Program-3D), and a three-dimensional code for isolated helicopter rotors in hover, 
ZAPR3D (Zonal Aerodynamics Program for Rotors-3D). The ZAP2D comparisons with large domain 
ARC2D solutions and with experimental data for a NACA 0012 airfoil have shown that the required 
domain size can be reduced to a few tenths of a percent chord for the low Mach and low angle of attack 
cases and to less than 2-5 chords for the high Mach and high angle of attack cases while maintaining 
solution accuracies to within a few percent. This represents CPU time reductions by a factor of 2-4 
compared with ARC2D. The current ZAP3D calculation for a rectangular plan-form wing of aspect ratio 
5 with an outer domain radius of about 1.2 chords represents a speed-up in CPU time over the ARC3D 
large domain calculation by about a factor of 2.5 while maintaining solution accuracies to within a few 
percent. A ZAPR3D simulation for a two-bladed rotor in hover with a reduced grid domain of about two 
chord lengths was able to capture the wake effects and compared accurately with the experimental pres- 
sure data. Further development is required in order to substantiate the promise of computational 
improvements due to the ZAPR3D coupling concept. 


i 
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NOMENCLATURE 


AR 

A-H 

Cf, C 4 , C m 


Cl> C d , C My 

Cq 

Ct 

j, k, f 

n 


R 


R, 

S 


V, v 
8 * 

<t> 

¥ 

A 

d 

y 


a 


Wing Aspect Ratio 
Grid domain labels 

Two-dimensional lift, drag and pitching moment 
Maximum two-dimensional lift coefficient 

Three-dimensional lift, drag and pitching moment 

Rotor torque coefficient 

Rotor thrust coefficient 

Grid indices 

Unit surface normal 

Local dimensionless rotor radius 

Domain boundary radius normalized by chord length 

Reynold’s Number 

Domain boundary surface 

Velocity field 

Rotor collective in degrees 

Surface velocity potential 

Azimuth angle 

Finite difference 

Partial derivative 

Surface vorticity 

Angle of attack in degrees 


Subscripts , Superscripts 

i, o Inner and outer boundary 

w Wake 

c Corrected 

g Geometric 
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1.0 INTRODUCTION 


The flow fields associated with full configuration, fixed- and rotary-wing aircraft are extremely 
complicated and remain a challenge for applied computational aerodynamics. While linear methods or 
potential flow models of nonlinear wake shape effects have become commonly used techniques for 
aircraft performance predictions and analyses, the nonlinear methods that incorporate more correct 
mathematical models of the actual flow physics— such as compressibility and viscosity— are being 
developed and applied by research scientists. Much progress has been made in solving the full poten- 
tial, Euler, and various forms of the approximate Navier-Stokes (NS) equations; however, because of 
the wide variation of physical scales of characteristic fluid phenomena from boundary layer thickness, 
shock thickness, or vortex core radius to wing span or fuselage length, a zonal method of analysis is 
generally projected for the entire aircraft flow field calculation. In such a zonal approach, the flow 
field is partitioned into zones which utilize approximate equations appropriate to the flow physics with- 
in each zone (e.g., Ref. 1 and 2). Although a zonal CFD code is more complex than, say, a NS code 
for the entire flow domain, required CPU time and computer memory should be significantly reduced 
and/or accuracy improved because of increased grid resolution where viscous or shock effects domi- 
nate. 


In the Phase I research effort, a new procedure for coupling potential and viscous flow calcula- 
tion schemes for a simple airfoil at low angle-of-attack was developed. In the Phase II work reported 
here, the two-dimensional concept is explored for the higher angle-of-attack range and is further devel- 
oped for transonic effects. The coupling procedure is expanded to three-dimensional wing simulations 
and then generalized for helicopter rotor simulations in the hover/climb flight condition. In the rotor 
procedures, rotational effects are added to the NS code and secondary blade and wake influences out- 
side of the NS domain are included in an outside iterative loop by a specialized potential flow hover 
code. The resultant two-zone methods, ZAP2D, ZAP3D and ZAPR3D (Zonal Aerodynamics 
Program-2D/3D and Rotor/3D), couples a potential flow panel method with a NASA Ames NS code, 
ARC2D and ARC3D respectively. 3 The coupling concept is based on the premise that any computa- 
tional method can only produce valid results within the approximations of the physical model employed 
in its construction. Therefore, the interfacing boundary surface between the potential and viscous flow 
code domains must be very nearly potential in order to achieve a successful (accurate) result. This is 
accomplished by utilizing an overlapped, velocity-coupling procedure with the unique feature that 
(beyond the first iteration) the potential flow surface panel strengths are obtained from the NS solution 
at a smooth inner fluid boundary. These fluid surface panel values are then used to compute the outer 
velocity boundary conditions for the NS calculation. The iteration between the potential flow solution 
and the NS solution continues in a closed loop until flow convergence, or allowable iteration limits, are 
obtained. 
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2.0 METHODOLOGY 


2.1 Two-Dimensional Investigation 

In this work, the coupling procedure is investigated and validated for two-dimensional high 
angle-of-attack and transonic airfoil simulations. The two-dimensional potential flow panel code, 
POT2D, has been fully coupled with the thin-layer, NS code, ARC2D. Both of these codes have 
already been extended to three-dimensional flows, VSAERO 4 and ARC3D, 3 respectively, and have 
been widely used and proven to be robust and computationally accurate flow simulation programs. 
The limits of the validity of the zonal procedure is established in the present work by comparing the 
ARC2D-alone simulation for large domain size with the fully-coupled ZAP2D simulation for small 
domain size for a range of angle of attack and Mach number. 

2.1.1 Zonal Concept 

In this zonal approach the flow field is partitioned into linear potential flow and viscous flow 
regions. The two computational zones overlap and the innermost coupling boundary is located in the 
computed flow field where the approximations inherent in both methods are valid; consequently, the 
interfacing boundary surface between the potential and viscous flow regions must be very nearly poten- 
tial in order to achieve an accurate result. 

The general coupling concept for the zonal method is illustrated in Figures 1 and 2. Figure 1 
shows the zonal representation of the physical flow field while the flow chart of the iteration scheme is 
shown in Figure 2. An inner boundary, S ; , which is measured by a radius, say Rj, is constructed to 
enclose a generally shaped body or airfoil and the fluid region near the body where viscous effects 
dominate. Although the body might be a complex shape, the inner boundary surface, S ; , would be a 
simple, smooth geometry. An additional outer surface, S„, which is measured by a radius, say R,,, also 
of simple shape, forms the computational domain of the viscous calculation. On the other hand, the 
potential flow domain extends from the inner boundary, S b to an infinite distance from the body sur- 
face. Consequently, the two computational domains overlap between the surfaces, S; and S 0 . Since S ; 
should be approximately a potential flow surface, the distance between Sj and S 0 should be minimized 
for optimum computational efficiency. Of course, the minimum overlap will be determined by numeri- 
cal error limitations. 

The iterative coupling of the zonal method proceeds as follows (see Figure 2). 

1. The calculation procedure is initialized by computing the potential flow associated with the 
body motion. An integral panel method is used for this computation. In this first step only 
the actual body surface is represented by discrete panels and the surface potential equation is 
solved directly by enforcing the Neumann boundary condition. 

2. The velocity field, V 0 , on the outer boundary, S 0 , is computed from the known potential flow 
surface singularity solution. 

3. This velocity field, V 0 , serves as the outer boundary values for the NS calculation. Once the 
NS solution for the flow inside S 0 is obtained, the velocity field due to the viscous solution, 
say Vi on S;, which is just outside the region of viscous effects, is also known. 
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4. This inner boundary velocity, Vj, can therefore be used to construct the known corresponding 
values of the potential flow panel strengths on S t ; consequently, all further potential flow 
calculations require only the paneling of the smooth inner boundary, S i; and not the actual 
body surface. In this way, numerical errors due to complex geometric modeling and calcula- 
tion of the outer boundary velocity field, V 0 , should be minimized and any contamination of 
the NS simulation due to this source of error should be avoided. 

5. This iteration loop, Steps 2-4, is repeated until the NS method is converged or iteration limits 
are achieved. Convergence should be measured by flow field and integrated force and mo- 
ment asymptotic behavior. 

2.1.2 Potential Flow Module 

As is well known, the potential describing an inviscid, irrotational, incompressible flow can be 
expressed by its boundary values. 5 Specifically, an integral solution of Laplace’s equation for the 
velocity potential can be derived by application of Green’s theorem. The resultant boundary integral 
equation for the surface velocity potential is composed of two integrals: 

1 . A body surface integral of a doublet kernel function of strength, </>, and a source kernel func- 
tion of strength, 3<£/3n, and 

2. A thin wake surface integral of a doublet kernel function of strength, A</> w . 

Of course, the wake strengths, A0 W , are related to the body surface potential, <f>, through the 
wake auxiliary Kutta condition and the body source strengths, d<j>ldn, are known by the tangent flow 
boundary condition on the body. Hence the integral equation for <t> can be solved directly as suggested 
by Morino. 6 In POT2D, the integral equation is discretized and transformed into an algebraic equation 
by representing the boundaries with flat panels and approximating the local panel source and doublet 
strengths with constant values just as in VSAERO. In the current code, the source strengths, d<j>ldn, 
are provided by the ARC2D normal velocity components on S-, and the matrix equation is solved for 
the corresponding potential doublet strength distribution, <j>, on S;. 

A wake model is used which implies a branch cut to infinity, thus canceling the vortex left at the 
upper and lower trailing edges of the S, boundary due to the truncation of the potential doublet panels 
(Figure 1). Additionally, an error may be introduced due to the truncation of the perturbation source 
strength at the rear of the inner surface. For this reason, the downstream outflow boundary remains at 
5 chords beyond the airfoil trailing edge for all grid domains tested. That is, the ARC2D solution near 
the outflow boundary at the coupling surface, S ; , is assumed to be near free stream conditions. 

Once the surface potential is computed, surface velocities are obtained by numerical differentia- 
tion and forces and moments are found by pressure integration. Further, off-body flow velocities are 
computed by integration of surface velocity source and doublet influence functions. The Karman-Tsien 
method 7 is used to correct the surface and field velocities for compressibility effects. 

2.1.3 Viscous Flow Module 

As noted earlier, ARC2D, originally developed by Pulliam and Steger 3 at NASA Ames, has been 
widely used in the computation of airfoil aerodynamics. The code has been continually improved by 
Pulliam and colleagues, 8- 9 • 10 and the details of the theory are well documented in these references. 
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The basic equations solved by ARC2D are the Reynolds-averaged NS equations written in strong con- 
servation-law form. These are simplified by including the standard thin-layer approximation for the 
viscous terms. The thin-layer approximation is known to be adequate for separation problems with 
subsonic or transonic free streams when the flow field associated with separation is convection-domi- 
nated. 

For turbulent flows, the well-known Baldwin-Lomax algebraic model 11 is used for turbulence 
closure. This turbulence model has been used in computing solutions for a wide variety of flow con- 
ditions and has been found to be acceptably accurate. 

The numerical algorithm is an implicit approximate factorization, finite-difference scheme, which 
can be either first- or second-order accurate in time. 

Over the years of development of ARC2D, unique features have been added to the code to im- 
prove its capability to accurately predict aerodynamic characteristics for airfoil calculations. One of 
these features is an implementation of a viscous-inviscid coupling. A point vortex correction model 
has been installed to reduce the solution dependency on the size of the outer boundary. Here, circula- 
tion on the far-field boundary is set based on the ARC2D computed C,. Using this velocity correction 
together with adjustments for compressibility, modified outer velocities are constructed. It is noted 
that this correction yields accurate results for domains larger than 4 chords in length. 12 

2.1.4 Boundary Conditions 

In ZAP2D, a parallel procedure is employed whereby the POT2D computed velocities with com- 
pressibility correction are used as far-field boundary conditions. At the outer boundary, velocities and 
pressure corrected for compressibility are specified and density is calculated from the pressure and free 
stream stagnation enthalpy. Since velocities, pressure and density are known, energy at the boundary 
can be determined. 

Since the POT2D formulation is based upon an incompressible flow assumption and the ARC2D 
formulation upon a compressible flow assumption, flow variables calculated at both boundaries must be 
compatible with each assumption. For this purpose, a simple method by Karman-Tsien, based upon 
the tangent-gas approximation, is used for the compressibility correction. At the inner boundary, 
velocity components calculated by ARC2D are converted to incompressible transpiration velocities for 
the POT2D representation. The pressure coefficient and velocity components calculated by POT2D at 
the outer boundary are in turn modified for compressibility for the ARC2D boundary conditions. 
These routines are alternated in the coupling procedure. 


2.2 Three-Dimensional Investigation 

The concept directly follows the previously described two-dimensional work, which employs an 
overlapping velocity coupling technique. As discussed, this method requires that the overlapping region 
is located where potential flow approximations are valid. Consequently, the actual wing surface and 
the viscous (or transonic) region near the wing are excluded from the overlapping boundaries. This 
zonal concept has been utilized to fully couple the widely used Navier-Stokes code, ARC3D, 3 with the 
potential flow panel code VSAERO. 4 
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The grid generation for this task was accomplished with the interactive code, HYGRID. devel- 
oped by Visual Computing, Inc. 13 This code was designed under the concept of Object Oriented Pro- 
gramming, and is capable of multi-block, hyperbolic grid generation about general three-dimensional 
configurations. 

The theoretical basis for the extension to three dimensions follows the two-dimensional version 
detailed in Section 2.1; therefore, further details are not required here. The validation of the resultant 
computer code, ZAP3D (Zonal Aerodynamics Program 3-D), will be discussed in Section 3.2. 


2.3 Three-Dimensional Wing/Vortex Investigation 

A wing/vortex interaction problem was investigated in order to test the requirement of capturing 
the vortex within the outer zonal boundary. The ability of the method to handle a vortex passing into 
the reduced zonal domain is fundamental to any rotor simulation. Two methods for this problem were 
investigated. One is to prescribe the known vortex core as fixed boundary conditions at the appropri- 
ate grid points within S 0 . This simple technique could eliminate the need for increased grid density 
near the vortex. The other is to use the S 0 boundary velocities perturbed by the vortex core. 

In either case, no changes were required to the basic potential flow module for the wing/vortex 
simulation. However, for validation purposes, a subroutine was added to ARC3D. This routine was 
utilized to test the accuracy of VSAERO to compute the induced velocities due to a vortex of known 
strength and location by the Biot-Sevart law. 


2.4 Hover/Climb Rotor Investigation 

The zonal method has been generalized for an isolated rotor in hover or climb. Since the experi- 
mental data by Caradonna and Tung 14 and calculations by Srinivasan, et. al. 15,16 and others 1718 were 
available for comparison purposes, a test case of a two-bladed, aspect-ratio-6 rotor with a NACA 0012 
section was chosen for this investigation. Some important issues were addressed, such as: 

a) Whether to modify ARC3D for a rotor flow simulation or use any readily available rotor 
code, 

b) Which type of grid system to be employed about the rotor configuration, and 

c) How to efficiently couple ARC3D, ZAP3D, and a rotor code to provide correct boundary 
conditions. 

Any viscous flow code alone may not properly calculate the flow of a lifting rotor unless its 
domain size is large enough to consider far boundaries as undisturbed ones, since the rotor induces 
significant velocities at large distances from the rotor. Therefore, specification of a no-flow condition 
at the inflow boundaries of a computational box, which is typically small for economy, poses a diffi- 
culty for the prediction method. This no-flow condition produces a closed box environment for the 
rotor in which the flow recirculates within the computational box. Here, a simple and economical 
computation program called HOVER, 19 which was developed by Summa using a vortex lattice method, 
is used to provide velocity boundary conditions at certain boundaries. (Nonlifting flow calculations can 
be done without using externally specified boundary conditions; namely using the freestream boundary 
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conditions.) For viscous flow calculations, various options were considered, which include the use of 
any rotor code available in the public domain and modification of ARC3D for rotational flow. Also, 
in the case of using ARC3D, the possibility of adopting either an unsteady approach for a rotating 
blade or a steady one for a fixed blade rotor with a rotating freestream was investigated. After some 
research, it was found that the modification of ARC3D for a blade-fixed representation could be 
achieved easily by adding rotational freestream terms in the subroutine RHS of ARC3D. 

2.4.1 Zonal Concept for HOVER 

For the three-dimensional rotor simulation, the zonal concept developed for wings in translation 
is generalized to include flow disturbances outside the NS domain. Since this general procedure is 
more complicated than described earlier, it is discussed in some detail here. The zonal representation 
of the hover/climb flow field is illustrated in Figure 3. Since this problem is axisymmetric, a single 
blade, say, the primary blade, is fully immersed within the NS domain, S„. While one might choose to 
simulate only the outer portion of the blade with the viscous code, it is hoped that the zonal coupling 
within S, will allow for the entire blade simulation and, hence, a root vortex. The effects of secondary 
blades and wake surfaces outside of S 0 are included by updating S 0 boundary conditions with a poten- 
tial flow method, the Analytical Methods, Inc. (AMI) HOVER code. 19 The technique of utilizing a 
potential flow hover calculation to initialize the finite difference domain boundary condition has been 
employed by other investigators. 14-20,21 However, the reduced outer boundary domain, S 0 , due to this 
potential/viscous coupling within S 0 , should allow for the closed-loop update of these outside potential 
flow influences. Such updating will be required to properly address viscous and compressible effects 
on hover/climb performance. 

The generalized zonal iteration scheme is shown in Figure 4. There are two iterative loops— an 
inside loop that uses the three-dimensional extension of the potential/viscous coupling procedure inves- 
tigated in Phase I and an outside loop to update the potential model of the flow disturbances outside of 
the NS domain which are driven by the viscous simulation. The iterative coupling of the general 
method proceeds as follows. 

1. The calculation procedure is initialized by computing a potential flow associated with the 
body motion. In fixed wing and body applications, an integral panel method developed by 
AMI, VSAERO, 4 would be used. For the rotor problems, programs HOVER and/or 
ROT AIR, 19 also developed at AMI, will be used for this step. 

2. The velocity field, say, K/ , on the outer boundary, S 0 , due to potential flow disturbances 
outside S 0 is computed. 

3. The velocity field, say, V' s , on the outer boundary, S„, due to potential flow disturbances 
inside S 0 is computed. 

4. The resultant velocity field, V 0 , serves as the outer boundary values for the NS calculation. 
Once the NS solution for the flow inside S„ is obtained, the disturbance velocity field due to 
viscous and compressible effects, say, V-, on S;, which is just outside the region of vis- 
cous/transonic effects, is also known. 

5. This inner boundary velocity, V ; , can therefore be used to generate known corresponding 
values of source strength, 3<j >/dn, and surface vorticity, y. These source and vorticity 
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strengths on S; are required for a new potential flow velocity update on S 0 , i.e., K s ' . Again, 
all further potential flow updates within S 0 require only the paneling of the smooth inner 
boundary, S ; . Therefore, numerical errors in Vj should be minimized and contamination of 
the NS simulation should be avoided. 

6. This inside S 0 iteration loop. Steps 3-5, is repeated until the NS method is converged or 
iteration limits are achieved. Convergence should be measured by flow field and integrated 
force and moment convergence. 

7. For the rotor case, the NS computed blade loading and wake development can then be used 
to update the potential flow simulation outside of S„. For example, secondary blades and 
wakes can be assigned new potential flow strengths and wake positions can be recomputed 
for fixed blade loading. 

8. With new blade and wake surface potential flow singularities obtained, an outside iterative 
loop is executed by returning to Step 2 and updating V £ . Once more, this outside loop can 

be continued until blade loading is unchanged or iteration limits are obtained. 

2.4.2 Potential Flow Module 

For the hover/climb case, the potential flow module included in ZAP3D for updating S 0 boundary 
conditions due to NS changes at S ; requires only minor changes. Specifically, the onset flow compo- 
nents as seen from the fixed-blade reference are extracted from the velocities at S, so that the Green’s 
function approach is applicable. Then, disturbance velocity updates to S„ due to the potential flow 
solution at S ; are added to the local onset flow components at S„. While these simple modifications 
take care of the secondary velocity coupling updates within ZAPR3D, the initial potential flow solution 
for the rotor flow field and also the secondary blade and wake influences outside of S; required the 
coupling with a full rotor simulation package, HOVER. 

HOVER is a "linearized" lifting-surface representation of the rotor airloads that is accomplished 
by a vortex lattice placed on the rotor planform area in the disk plane as illustrated in Figure 5. In 
HOVER, the influences of individual panels in the blade lattice are computed by quadrilateral vortex 
rings; therefore, the basic unknowns in the flow tangency equations are the panel ring vortex strengths, 
or, equivalently, panel doublet strengths. The program includes prescribed as well as relaxed wake 
calculations. If the user has supplied elastic blade properties, the elastic twist and bending deforma- 
tions and their impact on the rotor loads are computed during the program iterations. The user can 
also request a thrust coefficient and the program will adjust the collective setting through the prescribed 
wake iteration to obtain the thrust. The relaxed wake calculation then proceeds at fixed collective. 

Once a converged wake geometry is computed by the prescribed wake or relaxed wake options, 
inviscid forces and moments on the blade bound vortex segments are then evaluated in the usual way 
by applying the Kutta-Joukowski Law. Of course, the chordwise and radial pressure jump distributions 
are also calculated. 

Finally, with the sectional coefficient of lift distribution known from the lifting-surface calcula- 
tion, the profile drag and, hence, profile and total torque must be determined by falling back on empir- 
ical data. This reliance on empiricism is, of course, removed in the ZAPR3D representation, which 
includes the full viscous terms. 
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The discrete vortex filaments shed from the trailing edge of each blade represent the hovering 
rotor wake, which quickly separates into two parts for conventional rotors— an inboard sheet of weaker 
vorticity and an outer tip sheet that rapidly rolls up to form a very strong tip vortex. HOVER also 
includes a simple model of the tip vortex shedding across the blade chords shown in Figure 5, which 
improves the prediction of aerodynamic loading near the rotor tip. 

The overall wake structure in the HOVER program is illustrated in Figure 6, and consists of 
near-, intermediate-, and far-wake regions. The dimensionless axial coordinates at the start of the 
intermediate- and far-wake regions are ZFAR1 and ZFAR2, respectively. The near-wake region gen- 
erally includes four vortex passes below the generating blade and is the region of wake relaxation. 
The intermediate-wake region serves as a "buffer" zone between the near-wake filaments and far-wake 
model. In the far-wake, each helical vortex filament is continued as a semi-infinite cylindrical sheath 
of uniform vorticity. The far-wake velocity contribution is then computed by the equivalent source 
disc located at ZFAR2. 

Generally, the calculated thrust at the completion of the prescribed wake iteration for fixed blade 
collective is too high, and a relaxed wake iteration is required to compute realistic rotor performance. 
This is because the empirical equations are derived from experimental data of thrust, torque, and tip- 
vortex filament geometry. Consequently, prescribed wake programs demonstrate good correlation for 
integrated loads for a large number of conventional rotors since the prescribed wake constants are 
intimately related to the theoretical methods used to construct the empiricisms. However, these meth- 
ods are less successful when compared with known collective settings, which implies that although 
integrated performance is predicted successfully, local sectional loads may not be properly calculated. 
Certainly, if the wake itself deviates from the experimental data base for the prescribed wake constants, 
performance predictions could be significantly in error. 

The relaxed, force-free wake procedure includes the Scully code model 22 and the self-induced 
velocity contribution derived from the work reported by Widnall. 23 Additionally, velocity components 
along vortex filaments are computed in cylindrical polar coordinates. The radial, V R , and axial, V z , 
velocity components are integrated over a time step that is adjusted by the average tangential, V T , 
velocity components across the particular wake segment. Consequently, the final wake azimuthal 
gridding remains fixed throughout the relaxation iterations, and wake deformations are computed in 
these azimuthal planes. 

The calculated radial circulation distribution is analyzed at each prescribed or relaxed wake itera- 
tion to calculate the inboard extent of the tip vortex roll-up. During the relaxed wake iterations, if the 
radial position of the maximum circulation shed into the tip vortex changes by more than 4% of the 
radius, and the option to relax only the tip vortex geometry is active, then the wake segmentation is 
regenerated, growing or eliminating the necessary "inner sheet" filaments, based on prescribed wake 
constants computed from the current relaxed wake position. 

Finally, in order to utilize HOVER in the ZAPR3D procedure, a general velocity scan subroutine 
was developed that also interfaced with the required ARC3D boundary condition surfaces. In particu- 
lar, procedures for locating the wake penetration points through S; and also for summing the wake and 
secondary blade velocity influences outside of S ; were devised. Additionally, three-dimensional graph- 
ics displays were developed for visualizing and debugging the rotor velocity fields required in the 
ZAPR3D simulation. 
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2.4.3 Viscous Flow Module 


Major modifications in the viscous modules include those in ARC3D and the grid generation 
routines. 

For a rotor version of ARC3D, rotational ffeestream terms as a function of radius were added to 
the subroutine RHS. Because the grid system is constructed on a blade-fixed reference, no unsteady 
approach is necessary. Of course, if forward flight is of concern, then the grid should be generated at 
respective azimuthal angle positions to reflect the proper rotational velocities as well as advancing 
ones. 


Many different sets of boundary conditions, including a characteristic-type boundary condition 
procedure, were examined to find the effect of each boundary condition. All the boundary conditions 
are explicidy specified, in order to obtain the most consistent solution behavior. At the wall, a no-slip 
boundary condition is used. The density at the wall is determined by a zero th order extrapolation. The 
pressure along the body surface is calculated from the normal momentum relation. To ensure the 
continuity across the wake cut for the H-grid topology, the flow quantities are determined by averaging 
the flow variables from both sides of the singular plane. At die far-field boundary, velocities and 
density are extrapolated and pressure is assumed constant. At the upstream and downstream planes, 
velocity components are prescribed by HOVER calculations. Pressures at both planes are extrapolated 
and density is assumed constant. At the root boundary conditions, density and pressure are assumed to 
be same as these at the inner grid point. Contravariant velocity components are extrapolated. Finally, 
at all boundaries, the total energy is determined from the equation of the state. For the nonlifting case 
only, the upstream and downstream flow properties of the H-topology grid are swapped using a period- 
icity condition in the blade azimuthal direction. 

2.4.4 Grid Module 

A body-conforming, finite-difference grid of a cylindrical type has been used for one rotor blade 
(with a rounded-tip cap) because of the periodic nature of the hover flow field. In this grid setting, 
since no forward flight is concerned, the free stream rotates in the clockwise direction while the rotor 
is fixed at a given azimuthal angle. In contrast to the grid system (a C-C type) constructed for the 
three-dimensional, rectangular wing in translation, a C-H type grid system was chosen such that C- 
mesh planes along the radial direction from the rotor axis could better describe the rounded tip geome- 
try, and H-mesh planes along the streamwise direction could be easily duplicated. Once the H-type 
mesh was generated on the rotor, the grid planes at both leading edge and the trailing edge were rotat- 
ed to represent the upstream and the downstream of the rotor. Although smaller domains erected on a 
part of the rotor configuration may be used for zonal calculations, domains of about two chords away 
from the surface which include the entire rotor blade, were generated. On the other hand, a nonlifting 
flow calculation included a hemispherical domain which was generated for periodic boundary condi- 
tions to be applied at both inlet and exit planes. In general, the grid was clustered near the leading and 
trailing edges and near the tip region to resolve the tip vortex. It was also clustered in the normal 
direction with a geometric progression in order to resolve the viscous flow near the blade surface. In 
these cases, there are about 10 to 15 points in the boundary layer with a spacing of the first grid point 
from the surface equal to 0.00005 chord. This corresponds to a y+ value of 0(1) in order to capture 
the viscous sublayer effect. Also, the inboard plane near the axis of rotation was located at a radial 
station equal to one chord length and a full cosine distribution of grid points was employed with an 
equal distribution of grid points at the tip region. 
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A simple algebraic grid generation method was used. By using surface normals on the given 
surface points, prescribed far-field normals, and a geometric progression ratio which is calculated from 
the characteristic size of domain, the grid system of limited azimuthal range is generated. 
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3.0 COMPUTED RESULTS 


3.1 Two-Dimensional Validation 

The effectiveness of the coupling procedure installed in ZAP2D 24 is measured by demonstrating 
accurate and converged solutions for a range of angles of attack at subsonic Mach number in section 
3.1.1 25 and at angles of attack near the two-dimensional 2C^ point over a range of free stream Mach 

Numbers in Section 3. 1.2, 26 together with decreasing domain size. The impact of the coupling tech- 
nique on solution convergence characteristics and on required boundary size for a given accuracy in 
lift, drag and pitching moment has been investigated. In this work, the NACA 0012 airfoil was select- 
ed for the validation cases because of the experimental data base. 27 

3.1.1 Angle of Attack Validation 

The validity of ZAP2D effectiveness in angle-of-attack cases is demonstrated with the NACA 
0012 airfoil for decreasing domain size in comparison with the large-domain ARC2D simulation and 
with the experimental database. Calculations assuming a fully turbulent flow field were made with a 
Reynolds number of 3.0 million based on the chord length and the free stream velocity corresponding 
to a Mach number of 0.3. Geometric angles of attack from the experimental data by Harris, 27 ranging 
up to the stall angle, were specified in the computations with the corrections suggested by Harris. 

3. 1.1.1 Grid Domain Cases 

A total of six grid domains for the NS calculations were utilized to investigate the blockage ef- 
fects of the outer boundary on ZAP2D performance. Table 1 shows the grids and the corresponding 
case identifiers, where R 0 is normalized by chord length. A basic C-type grid around the NACA 0012 
airfoil was generated by an algebraic method where surface normal is maintained (Figure 7). The total 
number of grid points ranges from 12,025 (185 x 65) for the largest domain with R„=25 chord lengths 
(Grid A) to 6,549 (177 x 37) for the smallest domain with R„=0.12 chord lengths (Grid F). Succes- 
sive grids (B through F) were obtained by stripping off outer grid points from Grid A. The upper and 
lower airfoil surfaces each have 80 grid points for all grids, while the wake surface has 12 points for 
the largest domain. The grid generated is clustered near the surface with a normal spacing of 
0.000025 so that the effect of the viscous sublayer can be captured, and also clustered at both the 
leading and trailing edges with a spacing of 0.002. For the POT2D calculation, the location of the 
inner boundary surface, S, was held fixed at a distance approximately 10% of the chord from the 
airfoil surface. 


GRID 

MESH SIZE 

R. 

1/R„ 

A 

185 x 65 

25.0 

0.04 

B 

181 x 57 

5.0 

0.20 

C 

177 x 45 

0.5 

2.0 

D 

177 x 39 

0.17 

5.8 

E 

177 x 38 

0.14 

7.1 

F 

177 x 37 

0.12 

8.3 


Table 1. Summary of Grid Domains. 
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3. 1.1.2 Validation and Domain Size Effects 


ZAP2D and ARC2D (Version 195) were ported to the Stardent Titan and also to the IRIS 
4D/240GTX workstation. To investigate the effects of domain size (Grids A-F) on the integrated 
forces and moment computed by ZAP2D and ARC2D calculations for fixed input parameters were per- 
formed with Grid A as a baseline calculation for each angle-of-attack case. The full NS option was 
utilized and calculations were continued for 1,000 iterations. Originally, it was expected that post-stall 
behavior would be investigated with ARC2D and the effectiveness of the ZAP2D coupling could be 
measured for flow separation conditions. However, various attempts, including a time accurate ap- 
proach, smaller At, and various dissipation coefficients, have not yielded consistent results for such 
cases. This shortcoming of the NS simulation has also been recently documented by Raghavan et al. 2 * 
Consequently, the maximum angle-of-attack case reported here is a= 14.53°. 

ZAP2D converged solutions were obtained for all domain sizes and angles of attack up to 
a= 12.86°. At a= 14.53°, ZAP2D diverged for Grids E and F. While this behavior might be due to 
numerical error in the ZAP2D boundary condition updates or violation of the potential flow assumption 
at the inner boundary, ARC2D itself might also exhibit some instabilities for such small domains. 

The robustness and stability of ARC2D for small domains was checked by imposing boundary 
conditions computed by ARC2D with the large Grid A domain. In this study, it was found that accu- 
rate solutions within the small domains were obtained only if the turbulent viscosity was calculated up 
to the edge of the outer boundary where the magnitude of the vorticity was less than 1 % of the maxi- 
mum level around the airfoil. This is in contrast to the standard 75% of the domain size that is coded 
in ARC2D. With this change to ARC2D, the ARC2D computed loads for Grids C-E, subject to the 
outer and downstream boundary conditions provided by the Grid A ARC2D solution, exhibited devia- 
tions compared with the Grid A solution of less than 0.2% of C, and of about 5% of C d and C m . On 
the other hand, a solution could not be obtained for Grid F; hence, with the current coding in ARC2D 
which includes turbulent viscosity throughout the entire domain, the smallest domain size where a 
converged ZAP2D solution could be expected for a= 12.86° would be Grid E. In fact, for a= 12.86°, 
ZAP2D converged for all grids including Grid F; however, as noted above, for a= 14.53, ZAP2D 
diverged for Grid E. 

Figures 8(a) through 8(c) illustrate the effects of domain size on section lift, C ( , drag, C d , and 
moment, C m , coefficients, respectively. Each data line shows the calculation for a fixed angle of at- 
tack. Based on the work reported in Reference 35, ZAP2D is almost independent of domain size for 
outer boundaries larger than R„=5; consequently, ZAP2D calculations are illustrated in Figure 4 for 
Grids B-F (R„=5 to R o =0.12). The target values for C„ C d , and C m are then taken as the respective 
ZAP2D computed values for a domain of 5 chord lengths; i.e., for zero error, each trend line would 
exhibit no deviation from its initial value. 

As shown in Figure 8(a), the deviation of the ZAP2D computed C, is less than 1% from its 
target values for low angles of attack and has a maximum of approximately 4% for the high angle-of- 
attack cases. Again, at a= 14.53°, ZAP2D diverged for Grid E. 

As expected, the computed drag coefficient is more sensitive to domain size. In Figure 8(b), the 
deviation from its target values is less than approximately 30 drag counts for alphas less than 5°. The 
deviation for Grid E increases to a maximum of 75 drag counts at a= 12.86°. This numerical scatter 
has been traced to small errors in the computation of "near" field velocities adjacent to the outflow 
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boundary where potential flow panels are quite large. On the other hand, the deviation for Grid C at 
the higher alphas is caused by small velocity inaccuracies generated by ignoring the near-field compu- 
tation altogether or by a lack of accuracy in the potential flow solution. Regardless, the computed C d 
for the smallest is less than 25 drag counts in error from its target values for all angles of attack. 

The domain dependency of the computed C m , illustrated in Figure 8(c), is quite small (the scale is 
unusually expanded). For all angles of attack, the variation from target values is less than 0.005. 

Figures 9(a) through 9(c) compare the angle-of-attack trends of the ARC2D (Grid A) and ZAP2D 
(smallest domain) calculations with experimental data 36 for C t , C d and C m , respectively. The ZAP2D 
and ARC2D angle-of-attack trends are almost identical. Results for the two codes differ by less than 
3% in C,, less than 7 counts in C d , and less than 0.004 in C m . Once more, the scale of the pitching 
moment data in Figure 9(c) is unusually expanded for clarity. 

In general, the calculations are in good agreement with experiment, other than at a= 14.53°, 
where the experiment showed a stall phenomenon. It is possible that, by using Harris’ formula of 
a c =a g - 1.55 x C„ too much correction in geometric angle of attack was made in this nonlinear re- 
gion, which resulted in a smaller angle of attack than the actual stall angle. However, in Figure 9(a), 
the computed d(C f )/da is almost flat and the stall may occur just beyond this angle. (Additional calcu- 
lations revealed that the stall does occur between 15° and 15.5°.) Of course, the uncertainty of the 
turbulence modeling and the location of transition will influence the stall characteristics. Before the 
stall, the deviation of calculated C, from experiment ranges from 1.8% at Qf=-1 .81 0 to 8.3% at 
a= 12.86°. Evidently, the viscous calculation has not properly represented the actual loss in circula- 
tion due to the thickening boundary layer above approximately 12° angle of attack. In Figure 9(b), the 
computed C d is approximately 20 drag counts larger than the experimental values for alphas between 
4.97° and 12.86°. While approximately 30% of this error might be due to ignoring the laminar bo- 
undary layer at the airfoil nose, the remainder must be due to the underprediction of the boundary 
layer thickness effects. In Figure 9(c) the experimental C m trend with angle of attack is accurately 
predicted by ARC2D and ZAP2D up to the stall angle (note expanded scale). 

Finally, the drag polar is shown in Figure 10(a) and the pitching moment curve in Figure 10(b). 
Again, up to the experimental stall angle, the comparison of the computed and experimental values is 
reasonably good. It should be noted here that an assessment of experimental results for the NACA 
0012 airfoil also show significant data scatter 29 which might explain some of the differences in the 
comparison. 

Distributions of pressure coefficient computed with Grid F are shown in Figure 1 1 and compared 
with experimental data for a range of angles of attack. As can be seen, the calculations agree very 
well with experiment at moderate alphas. However, at the higher alphas, the computed pressures are 
too low compared with those of the experimental data. Again, apparently the boundary layer effects 
are not simulated exactly. However, all suction and pressure peaks are correctly predicted. Mach 
contours around the airfoil at a= 12.86° are shown in Figures 12(a) and (b). Figure 12(a) from the 
ARC2D-alone solution (Grid A) is almost identical to Figure 12(b) from the ZAP2D solution with Grid 
F. A substantial transonic and supersonic region exists near the leading edge (a maximum Mach num- 
ber of 1.2) in a region close to the body surface. Although some Mach numbers reach more than 0.5 
at the outer boundary, these are local phenomena and the compressibility correction used here is still 
considered reasonably valid. Figures 13(a) and (b) show pressure contours around the airfoil which 
agree very well. The normal pressure gradient near the leading edge on the upper surface is very large 
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and immediately downstream of the trailing edge the normal pressure gradient is negligible. This is 
true for the range of angles of attack investigated, and is confirmed by the streamline plot in Figure 14 
which is calculated using ZAP2D with Grid F at a= 12.86°. As shown, the streamline curvature ef- 
fects are important at high incidences near the leading edge and in the wake region where the flow 
tends to return to the free-stream direction. 

The numerical stability and convergence characteristics are shown in Figure 15(a) where the time 
histories of C, and C d are compared for 1,000 time steps for ARC2D (Grid A) and ZAP2D (Grid C). 
Beyond about 200 iterations, both solutions monotonically converge. The ZAP2D time history demon- 
strates the reduced magnitude of both C, and C d fluctuation levels due to a large reduction of the grid 
domain. Additional improvement in this behavior was obtained by incorporating an automatic schedule 
procedure in ZAP2D. In this scheme the first POT2D calculated boundary conditions are not updated 
until the rate of change of C, (AC,/At) in the ARC2D solution is less than 0.001. Once the boundary 
is updated by POT2D, the ARC2D calculation proceeds until AC f /At in the ARC2D solution is again 
less than 0.001. This procedure stops when both AQ/At is less than 0.001 and the L2 norm of the 
residual in the ARC2D solution is reduced by 3 orders of magnitude from the initial value. A total of 
four zonal updates results in a significant improvement in the solution history. In the range of 400 
iterations, the rate of change of C, and C d in the automated procedure is about half that of ARC2D and 
about the same as that of ARC2D at iteration 700. 

Finally, the time history of the L2 norm of residual is compared for ARC2D (Grid A) and 
ZAP2D (Grid C) in Figure 15(b). The ARC2D simulation shows a reduction of 3 orders of magni- 
tude, whereas the ZAP2D calculation converged to a reduction of about 4 orders of magnitude. It is 
natural that the ARC2D solution with Grid A should have a larger error of mass conservation than that 
of the ZAP2D solution with Grid C. Oscillatory behavior is observed in the L2 norm of residual in 
the ZAP2D calculation and is probably due to machine zero accuracy. Further, in the automatic 
scheduling, the density residual jumps abruptly at each update due to the sudden change of the bound- 
ary values. Although the calculation stopped after the density residual was reduced by 3.5 orders of 
magnitude, it is possible to allow more iterations to reach a lower residual. 

3. 1.1.3 Improved CPU Time 

For the low angle-of-attack case, a reduction of CPU time by a factor of 4 was reported in previ- 
ous work. 24 A reduction of grid points by a factor of about 2 and of iterations by the same factor were 
the result of using the automated zonal update schedule. For high angle-of-attack cases, a similar 
overall savings in CPU time was achieved. Total execution time on the Titan at a= 12.86° is 4,300 
seconds for ZAP2D with Grid C and 16,215 seconds for ARC2D with Grid A. The CPU times re- 
quired in these calculations were 0.00134 seconds/grid point/iteration for ARC2D and 5.5 seconds/ 
iteration for POT2D on the Titan; 0.00080 seconds/grid point/iteration were required for ARC2D and 

3.3 seconds/iteration for POT2D on an IRIS workstation. 

3.1.2 Mach Number Validation 

The objective of this work was to determine the range of applicability of this zonal concept for 
free stream Mach numbers through 0.84 at angles of attack near the two-dimensional C, point. In 
the following sections, the detailed analysis of the zonal simulation of the flow past an airfoil is com- 
pared with the NS alone simulation as well as with experimental data. 
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3. 1.2.1 Grid Domain Cases 


A total of eight grid domains have been utilized to study the sensitivity of the ZAP2D zonal 
concept to domain size with particular emphasis on compressibility effects. A description of the grids 
and their corresponding case identifiers are included below in Table 2. In this table and throughout 
this report, R 0 is the outer domain radius normalized by chord length. The total number of mesh 
points ranged from 12,025 (185 x 65) for the largest domain (outer boundary characteristic distance, 
R„ - 25 chord lengths— Grid A) to 6,919 (187 x 37) for the smallest domain (outer boundary character- 
istic distance, R„ - 0.1 chord length— Grid H). The upper and lower airfoil surfaces each had 80 grid 
points for all domains. A total of 12 mesh points was included in the wake region, for Grid A, be- 
tween the airfoil trailing edge and the downstream boundary at a distance of 25 chord lengths. To 
allow the use of the simple wake model in ZAP2D, the wake region for Grids B through H was mod- 
eled with 14 mesh points spaced from the airfoil trailing edge to a point 5 chord lengths downstream. 
This downstream boundary location remained fixed for all ZAP2D calculations to minimize the grid 
variation during the study. Also, the downstream grid spacing was held constant beyond approximate- 
ly one chord length from the airfoil trailing edge in an effort to minimize errors in the potential flow 
field velocity calculation. The ARC2D calculations were made exclusively on Grid A and were used 
as "target" values for the smaller domain, ZAP2D computations which were conducted on Grids B 
through H. The point vortex correction available in ARC2D 3 was not used in this study. This option 
has demonstrated a similar domain reduction capability in two-dimensional airfoil applications, but the 
ZAP2D concept has the advantage of possessing a direct three-dimensional extension. 


GRID 

MESH SIZE 

K 

1/Ro 

A 

185 x 65 

25.0 

0.04 

B 

187 x 57 

5.0 

0.2 

C 

187 x 54 

3.0 

0.33 

D 

187 x 52 

2.0 

0.5 

E 

187 x 48 

1.0 

1.0 

F 

187 x 45 

0.5 

2.0 

G 

187 x 40 

0.2 

5.0 

H 

187 x 37 

0.1 

10.0 


Table 2. Summary of Grid Domains. 


3. 1.2.2 Validation and Domain Size Effects 

Previous numerical studies to validate the basic ZAP2D zonal concept for an airfoil at a fixed, 
moderate angle of attack at low subsonic Mach number 24 and to extend the validation to a range of 
angles of attack up to stall were detailed in Section 3. 1.1. 25 Further validation of this zonal concept, 
included here, involves the flow about a NACA 0012 airfoil at Mach numbers 0.3, 0.5, 0.7, 0.8 and 
0.84. 26 The angles of attack for each Mach number (see Table 3) correspond approximately to the 
two-dimensional C t point as defined by the test data of Harris, 27 and as such represent a very chal- 
lenging set of conditions to test the coupling concept. In fact, it is assumed that these flight conditions 
will define the expected minimum benefit from this two-dimensional zonal method in terms of CPU 
time reduction. 
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Mach 

Number 

Apha 

(Geometric) 

Alpha 

(Corrected) 

mm 

14.86 

12.86 

wSmm 

9.86 

8.35 

0.70 

5.86 

4.75 

0.80 

2.86 

2.29 

0.84 

1.86 

1.67 


Table 3. Summary of Flight Conditions 


Since the inner fluid boundary surface, S ; , is to approximate a potential flow surface, care must 
be taken to ensure that the effects of the transonic flow conditions are minimized. Consequently, the 
local Mach number was monitored on the inner S, boundary at periodic intervals in the iteration se- 
quence. If the local Mach number achieved a level of 0.9 anywhere on the S, boundary, the code 
would automatically shift the indicia! location of Sj outward. For domains too small to satisfy this 
restriction, the inner S; boundary was located no closer than two grid contours inside S 0 . The choice 
of a "C" topology, of course, makes this type of adjustment very simple to implement. Further work 
would include utilizing an arbitrary shape for the inner Sj boundary. 

The effects of domain size (Grids A-H) on the computed integrated forces and moment have been 
thoroughly investigated for fixed ARC2D input parameters. Also, the ZAP2D results presented here 
include boundary condition updates at every iteration in the viscous/potential flow cycle. The solutions 
for Mach numbers 0.5 through 0.84 have been assumed converged when the solution residual achieved 
a third order reduction. The Mach 0.3 results exhibited the "stiffness” associated with low Mach 
number computations 41 and as a result were assumed converged when the solution residual achieved a 
fourth order reduction. 

Of course, in most cases a combination of force coefficient and residual reduction is a superior 
means of determining solution convergence. Here, it is only necessary to set some reasonable conver- 
gence criteria for comparing two numerical results. 

The lift, drag and pitching moment variation with domain size is presented in Figure 16. Again, 
all calculations presented are at the experimental C, values. The ARC2D "target” solution is the 

value located at an outer boundary radius of 25 chords (l/R o =0.04). The general trend for the lift 
variation with domain size (Figure 16(a)) appears to be consistent throughout the Mach number range 
tested. Some inconsistencies are present in the drag and pitching moment results in the region of Mach 
0.5 (Figures 16(b) and (c)). This behavior may be an indicator of the need for grid refinement for 
each flight condition to obtain optimum shock/boundary layer resolution. However, the validation of 
the zonal method by means of comparison with two related numerical algorithms is not compromised 
since the common NS component in ARC2D and ZAP2D is exposed to the same degree of local grid 
resolution. 
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The effect of compressibility at these high lift conditions upon the accuracy of the ZAP2D model 
is clearly shown as the freestream Mach number is varied from 0.3 to 0.84. To maintain the same 
level of accuracy (as measured by the difference between the ZAP2D integrated force and moment 
coefficients and the target values), the domain size must increase with increasing freestream Mach 
number. By requiring a deviation no larger than 3% (in absolute value) from the target value, only a 
20% chord domain is needed at Mach 0.3, while at Mach 0.84, a 5-chord domain is necessary (Figure 
17). Computing the lift to this accuracy resulted in a range in drag deviation from the target value of 
roughly 2% (10 drag counts) at Mach 0.84 to 10% (30 drag counts) at Mach 0.3. 

A plot of the experimentally determined maximum lift coefficient as a function of freestream 
Mach number, for the cases considered here, approximates the stall limit for the NACA 0012 airfoil 
(Figure 8(a)). To verify that the numerical results corresponding to the stall limit, further cases must 
be examined to determine the computed C, point. A 1.67° corrected angle of attack at Mach 0.84 

was the largest angle for which a converged (non-oscillatory) ARC2D result was obtained. Similarly, 
the work described in Section 3.1.1 had determined that the maximum corrected angle of attack at a 
Mach number of 0.3 is 12.86°. The computed maximum lift at the remaining Mach number points has 
yet to be verified. Presented in Figure 18(a) are lift results comparing ARC2D results on Grid A 
(R„=25) with ZAP2D results on Grids C and G (R c =3 and 0.2). As indicated, the ARC2D and 
ZAP2D (Grid C) results compare well with the data from experiment. The Grid G computations are 
included here to demonstrate the degradation in the small domain solution at the higher Mach numbers. 
The computed lift coefficient on Grid G at Mach 0.3 is within 2% of the ARC2D target value, while at 
Mach 0.5, the ZAP2D calculation is greater than 7% in error. It is fortuitous that the computed lift 
coefficient on Grid G at Mach 0.84 is nearly identical to the experimental value. 

Similar conclusions may be obtained when examining the drag and pitching moment variation 
with freestream Mach number at the flight conditions considered (see Figures 18(b) and (c)). The 
ARC2D and ZAP2D computed pitching moment coefficients are in close agreement with the data from 
experiment except at Mach 0.8 and 0.84 (Figure 18(c)). Here, an increased nose-down pitching mo- 
ment is apparent due to a computed aft shock position as compared with experiment. 

Further insight may be gained through examination of some critical flow features at a single 
flight condition. For this purpose a freestream Mach number of 0.7 and an angle of attack of 4.75° 
was arbitrarily selected. Contours of local Mach number are presented in Figure 19 for both ARC2D 
(Grid A) and ZAP2D (Grid C). The comparison between the two computed solutions is quite good for 
this challenging transonic flight condition. The ZAP2D result used for comparison corresponds to a 
3-chord outer domain radius (Grid C), which yielded a lift coefficient within 2.5% of the target 
ARC2D value. The drag coefficient was computed to within 2.8% (1 1 drag counts) of the target val- 
ue. 


The distribution of surface pressure coefficient, Cp, computed by ZAP2D for Grid C is com- 
pared with values from experiment in Figure 20. Deviations in surface pressure, between ARC2D 
(Grid A) and ZAP2D (Grid C), were not discernible to the scale shown in this figure. The correlation 
between ZAP2D and experiment is reasonably good over most of the airfoil with the exception of the 
shock location and downstream shock/boundary layer interaction region. This type of correlation is 
similar to that reported by Maksymiuk and Pulliam 12 for the NACA 0012 airfoil at Mach 0.8. Proper 
analysis of this transonic condition would require grid refinement which was beyond the scope of the 
present study. 
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The convergence characteristics are shown in Figure 21 where the time histories of C, are com- 
pared for ARC2D (Grid A) and ZAP2D (Grid C). Previous work at Mach 0.3 24 - 23 had indicated that 
the ZAP2D calculation was a stabilizing influence based on comparison of the convergence trends. As 
indicated here, the convergence history at Mach 0.7 is quite similar for the two codes. The previously 
observed convergent oscillatory behavior over the first 200 iterations has been replaced by highly 
damped oscillatory behavior over this time-step domain. 

The time histories of the density residual are compared for ARC2D (Grid A) and ZAP2D (Grid 
C) in Figure 22. The ZAP2D zonal concept demonstrates a faster reduction in residual compared with 
ARC2D alone. The ARC2D computation has achieved a third-order reduction in residual by iteration 
413, while ZAP2D requires only 281 iterations. 

3. 1.2.3 Improved CPU Time 

The initial work presenting the ZAP2D concept 24 had shown a saving in CPU time by a factor of 
4 for a low-Mach-number, two-dimensional airfoil computation. This represented a calculation on a 
grid domain radius of 0.14 chord lengths with an accuracy of 1% in lift when compared with the large 
domain result (Grid A). The work presenting the validation of this zonal concept at low Mach number 
for a range of angles of attack up to stall detailed in Section 3.1. 1 25 had shown a similar ratio. The 
results presented here indicate a saving in required CPU time by roughly a factor of 2 for Mach num- 
bers greater than about 0.5 at angles of attack near C, . In the present study, the ZAP2D zonal 
concept has demonstrated an increased rate of convergence over ARC2D alone and also a reduction in 
the grid domain radius from 25 to 3 chords (for Mach 0.8 or less), with less than 3% loss in accuracy. 


3.2 Three-Dimensional Fixed Wing Validation 

The objective of this work was to extend the previously documented two-dimensional zonal con- 
cept to three dimensions and to validate the technique for an aspect ratio 5 rectangular wing consisting 
of NACA 0012 airfoil sections. The following sections detail the validation process and include com- 
parisons between ZAP3D, ARC3D and available experimental data. 31 

3.2.1 Grid Domain Cases 

A total of six grid domains have been utilized to study the sensitivity of the ZAP3D concept to 
domain size for a rectangular wing at a single flight condition. A description of the grids is included 
in Table 4 below. 


GRID 

MESH SIZE 

R. 

1/Ro 

A 

100 x 39 x 23 

16.78 

0.06 

B 

100 x 39 x 20 

3.41 

0.29 

C 

100 x 39 x 19 

2.00 

0.50 

D 

100 x 39 x 18 

1.18 

0.85 

E 

100 x 39 x 17 

0.69 

1.45 

F 

100 x 39 x 16 

0.41 

2.44 


Table 4. Summary of Grid Domains. 
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The total number of mesh points ranged from 89,700 for the largest domain outer boundary 
characteristic distance (R„ - 16.78 chord lengths— Grid A) to 62,400 for the smallest domain (R„ - 0.41 
chord length— Grid F). The upper and lower airfoil surfaces each had 80 grid points for all domains, 
based upon a C-0 topology (see Figure 7). 

3.2.2 Validation and Domain Size Effects 

To validate the implementation of the ZAP3D coupling concept, a very detailed investigation of 
the induced velocity components on S ; , the resultant potential solution on S ; , and the induced velocity 
updates at S 0 are required. An artificial potential flow only zonal loop was constructed as utilized in 
the earlier two-dimensional work in order to validate the accuracy of the potential flow off-body veloc- 
ity calculation and the transpiration based coupling technique. By utilizing this method, special near- 
field (scan points that fall within a panel characteristic dimension) velocity parameters that are inherent 
in the VSAERO velocity scan procedure were sized to limit any induced velocity error to less than 1% 
for the expected S 0 and S, relative grid densities. 

As discussed earlier, since the "wake extension" of the S; boundary includes no source terms, the 
outflow boundary should theoretically be located far enough downstream such that ARC3D obtains 
nearly ffeestream conditions on S ; at this location. However, for outflow boundaries located far away, 
the wake will eventually impact S, or even S c for onset flows at non-zero angles of attack. A compari- 
son of the velocity components on S, as computed by ARC3D and the standard ZAP3D simulation 
(Figures 23(a) and (b)) indicates the effect of wake impingement on the coupling surfaces. This condi- 
tion violates the potential flow approximation of the ZAP3D zonal concept at S r Consequently, the 
computed ZAP3D wake behavior is completely incorrect. That is, ARC3D information involving the 
wake deficit is passed to the Sj boundary in the form of a normal velocity boundary condition only, and 
there is no mechanism in the current model to accurately transmit any vorticity delay information via a 
potential flow surface to the outer boundary S 0 . 

Furthermore, because the grid density is not maintained to capture the wake correctly for the 
entire distance downstream, numerical diffusion will certainly lead to errors in the Navier-Stokes ve- 
locities on the Sj boundary. While not thoroughly investigated here, such local errors apparently may 
not significantly affect the computed wing forces and moments in the ARC3D-alone simulations, but 
can drive the ZAP3D simulation to a different converged result because of the potential flow update of 
the boundary conditions at S 0 . 

One technique that can be utilized to minimize the wake-impingement problem is to align the 
downstream grid with the ffeestream direction. This, however, would require rotation and redistribu- 
tion of the grid for each angle of attack. Clearly, other more practical corrective measures that do not 
involve such geometric changes are preferred. 

The technique adopted here, which avoids modification of the volume grid, is to shorten the 
downstream extent of the effective S, boundary while maintaining the entire S„ boundary. Of course, 
as S; is shortened, the "wake extension model" which is required in ZAP3D is also brought forward 
within the S 0 domain. This modification avoids the problems associated with the wing wake impinge- 
ment on S;. It does, however, imply some additional error of unknown but probably small magnitude, 
since the "wake extension" of S ; does not presently include any source term that should be required to 
continue ARC3D disturbances from ffeestream conditions at the S ; terminus. A comparison of the 
ARC3D and ZAP3D velocity components with the shortened Sj surface (Figures 24(a) and (b)) indi- 
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cates a substantial improvement in the streamwise velocity component. The vertical component shows 
some improvement over the wake region, but is slightly worse near the airfoil trailing edge station 
(note the scale differences between Vx/a and Vz/a). For the present case, evidently the potential errors 
due to the neglect of any source terms on the wake extension appear to be small; therefore, this meth- 
od has been incorporated in all the calculations presented below. 

As previously mentioned, domain size effects for ZAP3D and ARC3D were obtained over a 
range of R 0 values from 16.78 to 0.441 chords (see Table 4). A single flight condition of a = 8.75°, 
Mach 0.12 and R, 1.2 million was examined for an aspect ratio 5 rectangular wing. 3 ' The ZAP3D 
simulation consisted of potential flow outer boundary velocity updates at every 250 ARC3D iterations. 
The S; boundary was located a distance of 0.14 chord from the wing for all ZAP3D calculations 
(/= 14). Furthermore, based on the previously described wake impingement study, the inner S; bound- 
ary was shortened to a location just beyond the wing trailing edge (x/c=1.02, j = 85). All of the 
ZAP3D and ARC3D calculations achieved a residual reduction of 3-4 orders. 

Figure 25(a) illustrates the ZAP3D and ARC3D lift variation with domain size. As indicated 
by the greatly expanded vertical scale in this figure, the ZAP3D three-dimensional lift C L for R o =0.69 
(1/R„=1.45, Grid E) is less than 3% higher than the target ARC3D C L with R 0 = 17. In this case, the 
shortened S, in ZAP3D is very effective in maintaining the ARC3D target C L for reduced domain size. 
Still, based on the two-dimensional results, 23 which indicated less than 2% error in Cl for R„ values as 
low as 0.14 chords, smaller domain sizes than achieved here should be possible. It may be that the 
three-dimensional wake impingement on S„ that still occurs in these calculations is limiting the amount 
of domain reduction that is possible for a given error limit. Further studies of the effects of the down- 
stream location of the outflow boundary and the grid density in the wake region should be investigated 
as a possible source of error that was not present in the two-dimensional simulations. Figures 25(b) 
and (c) illustrate the computed drag coefficient (C D ) and pitching and moment coefficient (C My ) varia- 
tion with domain size. As illustrated, for ZAP3D domain size reductions down to 1 chord length, 
there are fewer than 10 counts variation in C D and less than 0.0015 change in C My . 

Because it exhibits the least error of the smaller domains in the computed C L compared with the 
ARC3D target value (less than 0.5% error), the ZAP3D result at an outer radius of 1.18 chords (Grid 
D) was selected here for further detailed analysis. A comparison of the symmetry plane chordwise 
pressures between the large domain ARC3D result (R„ = 17) and the ZAP3D calculation is presented in 
Figure 26. There is no discernable difference to plotting accuracy in the two simulations. The com- 
puted spanwise loading distributions by ARC3D and ZAP3D are compared with the wind tunnel data 
in Figure 27. Again, the ARC3D and ZAP3D results are nearly identical and compare very well with 
the experimental data. The history of solution residual for the ARC3D large domain and ZAP3D small 
domain is presented in Figure 28. The spikes in residual for the ZAP3D calculation are associated 
with the outer boundary velocity update (potential flow update every 250 ARC3D iterations). The 
relative improvement in residual reduction over ARC3D is due to the reduced distance over which 
information must travel. These phenomena were also noted in the two-dimensional work, 24 Section 
3.1. 


The history of total lift coefficient is shown in Figure 29. The improvement in C L convergence 
over the first few hundred iterations that was noted in the two-dimensional work is not apparent here. 
Additional ARC3D large domain data points are required to be conclusive, since there is lift coefficient 
information only at every 250 Ih iteration. 

As noted earlier, the final deviation in lift coefficient between ZAP3D (R„=1.18 Grid D) and 
ARC3D (R„=17 Grid A) is less than half a percent. To demonstrate the importance of multiple vis- 
cid/inviscid iterations, an ARC3D calculation that is initialized with the potential flow solution is 
shown in the inset in Figure 29. As illustrated, the ARC3D lift coefficient, given a potential flow 
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initial condition, deviates from the multiple update ZAP3D calculation (and ARC3D large domain) by 
about 3%. 

Finally, a comparison between the ZAP3D computed pressure coefficient and the values from 
experiment 31 are presented in Figures 30(a), (b), (c) at inboard, mid, and outboard wing span locations. 
Overall, the comparisons in the pressure distributions are quite good. 

3.2.3 Improved CPU Time 

A general statement regarding the saving in CPU time by the three-dimensional zonal method is 
premature at this time. The outflow boundary location and downstream grid density effects on the 
reduced domain calculations need to be further investigated to achieve the minimum grid domain sizes 
that should be possible. Furthermore, higher angles of attack and higher Mach numbers need to be 
studied before the issue of CPU time is resolved. 

Nevertheless, based upon the convergence behavior shown in Figures 28 and 29, the current 
ZAP3D calculation at an outer domain radius of 1.18 chords represents a speed-up in CPU time over 
the ARC3D large domain calculation by a factor of about 2.5. This improvement is obtained for less 
than a 0.5% deviation in C t , 10 counts change in C D , and 0.0015 variation in C^. Since the two- 
dimensional studies achieved a speed-up by a factor of 4, it is expected that further reductions in the 
required computational domain size can be achieved with sizable improvements in CPU time. 


3.3 Three-Dimensional Wing/Vortex Validation 

To validate the scheme for vortex-induced lift on a wing, a simple numerical calculation was 
made. A vortex -generating wing (aspect ratio 5) with a NACA 0012 airfoil section was mounted 
vertically with 8° of yaw and located upstream of a horizontally-mounted wing. This downstream 
wing was geometrically identical to the upstream one, but was set at zero angle of attack. Other flow 
conditions were M=0.12 and Re=1.5 million, based upon the chord length. The flow on the investi- 
gated wing was assumed to be fully turbulent and the Baldwin-Lomax algebraic turbulence model was 
used for turbulence closure. Figure 31 shows a typical location and shape of the far boundary sur- 
rounding the wings (R„= 18). 

As explained in Section 2.3, the upstream wing and shed tip vortex will modify the boundary 
conditions at S 0 . Various outer domain boundaries were investigated in order to test the ability to 
capture these effects within ARC3D. After some study of boundary locations, five different domain 
sizes were thoroughly investigated, with boundaries ranging from a large radius outer boundary which 
completely encloses the wings (R„=5.8) to a small radius outer boundary just outside the investigated 
wing (R o =0.08). Up to about 90,000 grid points were used for the largest domain with 98 points 
along the streamwise direction, 43 points along the spanwise direction, and 21 points in the normal 
direction. Smaller domains were achieved by removing outer layers of the C-type grid. In all domain 
sizes studied, the downstream exit boundary remained at a fixed location. At this exit boundary, con- 
ditions are internally computed from ARC3D. 

Figure 32(a) illustrates the computed spanwise lift coefficient of the downstream wing for the re- 
spective outer boundaries. The computed span loadings illustrate the effect of the upwash and down- 
wash induced by the tip vortex located at Y = 2.0 (80% of span length) and Z=0.5. Among five 
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ARC3D/VSAER0 solutions, only two— with the outer boundary at R o =0.4 and 0.24— are regarded as 
valid because only these two show reasonable comparisons with the VSAERO-alone solutions. The 
solution with R„=0.08 apparently has the boundary too close to the viscous region around the investi- 
gated wing. The cases with R o =0.68 and 5.8 contain the tip vortex within the domains, but the vortex 
effects were not properly accounted for in the calculations because of numerical diffusion due to the 
course grid density for these larger domains. Clearly, if the vortex is present within the calculation 
domain, the grid should be dense enough to capture the vortex from the outer boundary condition and 
maintain its effect without dissipation, or it should be modeled as a potential prescribed vortex filament 
which is carried throughout the domain. 

In another study to thoroughly resolve the differences obtained with ZAP3D and ARC3D initial- 
ized by VSAERO for the vortex/wing interaction test case, the velocity updates at the S 0 boundary 
were replaced with an exact line vortex calculation. In this way, any numerical errors arising in the 
ZAP3D/VSAERO updates associated with the boundary velocities would be eliminated, and the 
ARC3D simulation could then be investigated for grid convergence. Again, the NACA 0012 wing was 
utilized with a vortex of known strength passing above the wing. The wing had a chord of 1 and a 
semispan of 2.5, and the infinite line vortex was located at Y v , Z v = (2.0, 1.5). The vortex strength 
was arbitrarily set to T=0.5, and symmetry conditions were enforced on the X-Z plane. Further, the 
wing angle of attack was set to zero so that any lift was entirely induced by the influence of the line 
vortex. 

Figure 32(b) illustrates the computed spanwise lift coefficient for several grid cases. The VS- 
AERO calculation is also shown for comparison and is considered the target solution, because a poten- 
tial flow approximation should accurately model this particular case. Although the figure includes 
calculations for two different domain sizes (R o =1.02 and R 0 =1.7), grid density studies have been 
completed only for the R 0 = 1.7 boundary. In this case, the vortex penetrates the NS domain between 
Si and S„. Furthermore, for all calculations, the vortex influence is included as velocity boundary 
conditions at S„ only. With R 0 = 1.7, three different grid densities were used to calculate the flow (98 
x 39 x 19, 98 x 39 x 31, and 98 x 38 x 51). In the denser grid cases, grid lines were added between 
the f = 19 and f = 18 surfaces such that the tip vortex passing between two surfaces could be more 
accurately resolved. The comparison with VSAERO shows that in the original grid (98 x 39 x 19) sec- 
tional lift was overpredicted on the inboard wing and underpredicted on the outboard wing. In this 
case, the vortex effects are incorrectly dissipated within the viscous solution. In the latter two grids, 
there is less than a 3% variation in the spanwise lift. Because the spanwise lift also compares very 
well with VSAERO for these two cases, the viscous effects are believed to be correctly captured. 
Hence, if the grid density near the vortex is sufficient, ARC3D can successfully simulate the 
vortex/wing interaction where the vortex penetrates S„, and the vortex effect is provided by velocity 
boundary conditions on S 0 only. For the rotor cases, this ability is a requirement in order to utilize the 
ZAP3D concept. 


3.4 Hover Rotor Validation 

3.4.1 A Nonlifting Rotor Case 

The nonlifting flow may be considered self contained when the flow field at the outer boundary, 
S 0 , lies undisturbed. In order to achieve this state, the outer domain, S„, should be located at least 10 
times the chord length beyond the rotor. Then, it is sufficient to impose the rotational freestream 
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boundary condition at the S„ boundary without using any externally provided boundary conditions; 
thus, NS alone calculations are possible. Because of the axisymetric flow, periodic boundary condi- 
tions (see Section 2.4.3), are imposed and the calculations are performed for only one blade. The grid 
at the inlet and exit planes maintains the exact symmetry and no interpolation was required to swap the 
flow variables in these planes. 

As a preliminary case to verify the rotor simulation, calculation of a nonlifting, two-bladed, aspect- 
ratio-6 rotor with a NACA 0012 section was carried out. Experimental data 14 and other NS calcula- 
tions 15 ’ 1 * are available for a tip Mach number of 0.52, and a Reynolds number (based on the tip speed 
and chord) of 2.32 million. At this Reynolds number, the flow region on the surface can be regarded 
as turbulent over the entire blade. For a viscous calculation, a domain size of R,=20 with about 
125,000 grid points of 91 points along the streamwise direction (40 points on the blade), 39 points 
along the spanwise direction, and 35 points along the normal direction were used. A partial view of 
the grid generated is shown in Figure 33. 

A typical ARC3D solution for the Cray-2 required about 2,000 time steps to reduce the residuals 
by 0(-3) with a CPU time per time step per grid point of 35 ms. The time step of 1.00 was used with 
the Jacobian correction to accelerate the convergence to steady state. The maximum CFL number cor- 
responding to the step size was 1.7. The calculated pressure coefficients, based upon the sectional 
dynamic pressure, are shown in Figures 34(a) through 34(e), and compared with the experiment. 
Although relatively coarse grid was used— compared with those in other similar calculations— results 
are in very good agreement. Consequently, the modifications to ARC3D for the rotor simulation are 
correct for this case, and the more interesting lifting case could be investigated with the ZAP3D zonal 
coupling. 

3.4.2 A Lifting Rotor Case 

The HOVER program installation on Unix machines and its modifications were verified by com- 
paring the wake geometry and integrated performance for an untwisted, rectangular planform, two- 
bladed model rotor of aspect ratio 13.7, which was reported in Reference 14. The measured wake 
geometry data was reported to be in agreement with the prescribed wake equations of References 18 
and 19. Hence, the calculated wake geometry should remain relatively fixed during the relaxation 
iterations in the HOVER module. 

In the HOVER initial solution, the full wake was relaxed and calculations were started with 
Landgrebe’s prescribed wake geometry. The computed collective was adjusted to obtain a thrust coef- 
ficient equal to that measured in the test. The computed wake geometry converged to a unique posi- 
tion regardless of the starting geometry. 

Figures 35(a) and (b) compare the measured and calculated tip vortex geometry for five relaxed 
wake iterations for C T =0.0037. The calculated geometry is well within the scatter of the experimental 
data, 11 and there is very little movement of the tip vortex through the relaxation iteration. Figure 36 
illustrates several constant azimuthal cuts through the inner sheet and tip vortex for the last relaxed 
wake iteration. For ^=90° and 180°, there are no differences in these geometries when compared 
with the second relaxation iteration. Moreover, the inner sheet does remain linear as expected with 
only a slight roll-up indicated at the outboard end where the vorticity is extremely weak. 
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The measured trend of thrust coefficient, C T , with collective, 0 75 , and torque coefficient, C Q , is 
also predicted quite well as shown in Table 5. 


MEASURED 

PREDICTED 

CD 

tn 

C T 

C Q x 10 3 

CD 

^4 

cn 

■ C T 

V* iq3 

5 ± h° 

.0018 

.109 

5.2° 

.0018 

.117 

8 ± h° 

.0037 

.253 

8.5° 

.0037 

.250 

12 ± 5° 

.0056 

.493 

12.1° 

.0056 

.439 


Table 5. Detailed Comparison of Measured and Calculated Hover 
Performance for the Ames Untwisted Rotor. 


The difference in the C Q for the 12° case is probably due to separations that occur at the lower 
test Reynolds number. The corresponding performance map comparison is shown in Figure 37 where 
the usual +2% experimental error band on C x for fixed C Q has been applied. These results verify the 
accuracy of the Unix installation and modifications of the HOVER code. 

As noted earlier, the coupling of the HOVER module within ZAPR3D is essential in order to 
complete the rotor simulations. In ZAPR3D, the HOVER module must first provide the converged 
rotor blade loads and associated rotor wake vorticity and geometry. Again, since the vortex-lattice 
method inherently includes an infinite computational domain, this computed rotor flow field should 
facilitate large reductions in the required viscous computational domain. The HOVER module then 
passes boundary condition information to the ARC3D rotor module by velocity scans at appropriate 
off-blade boundaries. Specifically, the entire HOVER flow field is utilized to compute the initial S 0 
boundary conditions for ARC3D while the HOVER flow field outside of Si is utilized to provide a part 
of the S„ boundary conditions for continued ZAPR3D iterations. During the ZAPR3D iterations, the 
completion of the boundary condition updates at S 0 are accomplished by velocity scans due to the ZAP 
potential flow representation of the inner surface coupling boundary, S;. 

In order to complete this work, the procedure for sharing S-, and S„ boundary information was 
generalized for the C-H grid topology used for the rotor simulation. Whereas the S; boundary for the 
wing cases were transmitted by specifying a single constant-f surface, the rotor case require transmit- 
ting several restricted grid surfaces. Additionally, a general procedure was also designed and coded 
to identity the HOVER wake intersections with S ; and the portions of the wake contained within the S ; 
volume. With the geometric work complete, a generalized velocity scan procedure was written for the 
HOVER module that would automatically compute both of the required ZAPR3D velocity updates at 
S„. Figure 38 compares the HOVER velocity scans at a ZAPR3D boundary located just upstream of 
the primary blade for the rotor case that is described below. The downwash predicted by HOVER for 
the entire flow field is illustrated in Figure 38(a), while Figure 38(b) illustrates the downwash for all 
influences outside of the S, volume. As expected, there is a reduction in the downwash due to the 
elimination of the wake effects within S s . While no comparisons with experimental data are possible 
for such a velocity decomposition, detailed analysis of several rotor calculations have been utilized to 
confirm the HOVER velocity scan routine. 


24 




With the HOVER velocity scan module complete, lifting rotor simulations by the modified 
ARC3D program and ZAPR3D could be researched. The test case for these viscous calculations is the 
same rotor that was investigated in the nonlifting ARC3D simulations. Numerical calculations for the 
hover condition are presented for a single collective setting of 8 degrees, a tip Mach number of 0.44, 
and a Reynolds number (based on tip speed) of 1.92 million. These conditions correspond to the 
available experimental test conditions. Again, in all calculations, the flow on the surface was assumed 
to be fully turbulent and the Baldwin-Lomax algebraic model was used for turbulence closure. 

While ARC3D only calculations cannot be obtained for the lifting rotor case, a preliminary 
ARC3D calculation that utilized periodic boundary conditions at the inlet and exit planes and HOVER 
provided boundary conditions at the far boundary located below the rotor was investigated in order to 
provide a baseline for the ZAPR3D simulations. For this purpose, a grid of an ellipsoidal shape (mi- 
nor axis along the rotational axis) was generated with 91 points in the rotational stream direction, 39 
points in the radial direction, and 35 in the normal direction to the blade surface. As usual in the work 
reported here, the normal grid spacing was clustered near the blade in order to achieve y+ values of 
0(1) for the grid points just off the surface. Because this rotor has two blades, the grid includes a total 
azimuthal angle sweep of 180 degrees. Also, in order to avoid interpolations of flow quantities at the 
periodic planes, wakecut planes at upstream and downstream locations with respect to the blade surface 
grid, which is itself at 8 degrees incidence, were gradually elevated and lowered respectively to match 
grid points at inlet and exit boundaries. Of course, this grid generation scheme is only required for 
calculations using periodic boundary conditions. For general boundary conditions which are employed 
for other lifting calculations illustrated here, please refer to Section 2.4.3. 

This large domain ARC3D calculation with periodic/HOVER boundary conditions achieved 0(3) 
residual reduction in 2000 iterations. While convergence was obtained, the grid was apparently too 
sparse throughout the volume to accurately resolve the entire flowfield wake effects, since the calculat- 
ed blade pressure distribution did not reflect the wake induced downwash or blade passage effects cor- 
rectly. Figure 39 illustrates a typical chordwise pressure distribution (Cp’s are based on the sectional 
onset flow dynamic pressure) comparison with the experimental data (and another calculation discussed 
later). As shown in this figure, the periodic/HOVER boundary condition solution does not compare 
very well with the data on the upper surface and also behaves poorly at the trailing edge. 

Rather than pursue vastly increased grid density that is required to resolve the wakes for the full 
180 degree domain, a reduced zone grid (Zonl) was generated that encloses the primary blade with 
upstream and downstream boundaries located approximately 1 chord (R„= 1) from the leading and 
trailing edges respectively. The grid included a mesh density of 89x39x35 for a total of about 121,000 
points. A view of several selected planes for the Zonl grid is shown in Figure 40. In this case, 
ZAPR3D is used for the simulation but no inside potential flow updates from Si are requested (i.e., 
ZAPR3D (ITI=0)). In this way, the HOVER solution provides the velocity conditions at all bound- 
aries for a single ARC3D converged simulation. The computed chordwise pressure coefficient at the 
68% radius is illustrated in Figure 39. Comparison with the experimental data is very good at this 
location and is clearly much improved over the periodic boundary condition case. The trailing edge 
behavior of this computation is also improved over the earlier result; however, comparisons with the 
experimental data at other locations indicate that the wake effects, while improved, are still not com- 
pletely resolved. This is illustrated in Figure 41 which compares the ZAPR3D (ITI=0) and HOVER 
computed radial distributions of the section thrust coefficient. Since the HOVER simulation compares 
accurately with the experimental data (CT=0.0048), the passing tip vortex effects have clearly not 
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been captured correctly in the viscous calculation. It also appears that conditions near the root may 
need to be improved as well. 

In order to capture the tip vortex correctly, a new grid (Zon2) was generated, and the full 
ZAPR3D simulation utilized. A perspective view of the new zone boundaries for the coupled calcula- 
tion is shown in Figure 42 while a quantitative display of the relative positions of the inner and outer 
boundary surfaces are shown in Figure 43. In this case, the outer boundary is located approximately 2 
chord lengths from the blade (R>=2) and the inner boundary approximately I chord length. For this 
new grid, Zon2, the grid density has increased to 89x79x55 for a total of about 387,000 points, or just 
over three times that used in Zonl. The radial density in Zon2 on the inboard half is the same as that 
used in Zonl but the outboard half density is much finer. Further, since the tip vortex trajectory was 
already known from the HOVER calculation, the grid density in Zon2 was increased in this region. A 
smaller geometric expanding progression was also used in the normal direction for Zon2. Finally, the 
wake cut distances at both the root and tip were maintained in the hope of improving the solution near 
the root. 

Two ZAPR3D simulations have been undertaken for this fine grid domain. First, a calculation 
that includes no inside loop updates, ZAPR3D (ITI = 0), was completed and required 2500 viscous 
iterations to achieve 3 orders reduction in residual. Second, a calculation that includes an inner surface 
potential flow update every 250 iterations, ZAPR3D (TTI=7), was started and continued for approxi- 
mately 1700 viscous iterations before running out of the requested computer time. Unfortunately, this 
calculation had achieved a residual reduction of only 2 orders at this time and because of time con- 
straints has not been continued at this time. Figure 44 illustrates the convergence history for these two 
cases for the first 2000 viscous iterations. Comparison of these characteristics with those obtained in 
the 2D and 3D wing calculations indicate that while the inner surface zonal update is being correctly 
modeled, the domain size needs to be further reduced in order to demonstrate the convergence trend 
that should be expected. Since the rotor case requires the relative comparison between ZAPR3D solu- 
tions rather than the comparison between small domain ZAP3D and large domain ARC3D solutions 
that were possible in the fixed wing simulations, the current trend reported in Figure 44 should be 
expected until the domain size is much smaller. 

Figure 45 compares the HOVER, ZAPR3D (ITI=0), and ZAPR3D (ITI=7) predictions of the 
radial thrust loading. Clearly, the ZAPR3D (ITI=7) solution is not yet converged as noted above and 
will undoubtably match the ZAPR3D (ITI=0) solution when finished. Furthermore, the ZAPR3D 
results for the increased grid density has more correctly captured the rotor wake effects than those 
reported in Figure 41. The ZAPR3D (ITI=0) maximum sectional thrust value does compare favorably 
with the HOVER value but the peak value is located further outboard in the viscous case. Additional- 
ly, the ZAPR3D loading near the root indicates that the boundary conditions imposed there may not 
yet be consistent with the HOVER simulation. These differences and the radial slope change in the 
ZAPR3D thrust loading at the mid-radius location may be due to abrupt grid distribution changes. 
More detailed grid studies will need to be investigated in order to resolve these issues. 

Finally, Figure 46 compares the converged ZAPR3D (ITI=0) computed pressure distribution for 
this rotor with the experimental data. Again, it is expected that for this domain size, the completed 
inside loop simulation, ZAPR3D (TTI = 7), will produce essentially the same computed pressure distri- 
butions illustrated here. Clearly, the wake effects are more correctly captured in this calculation than 
the earlier simulations (compare with Figure 39). In fact, this calculation compares very favorably 
with the experimental data with the exception of the suction peak near the root and the pressure on the 
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suction side near the trailing edge at the 98 % radius location. This later trend is most likely due to 
the tip edge vortex that will probably require additional grid refinement. 

Future ZAPR3D work should concentrate on detailed grid refinement studies and systematic 
reductions in the outer boundary domain size. Undoubtably, such studies would remove some of the 
limitations of the simulations obtained so far and would also facilitate the detailed validation of the 
effectiveness of the inner boundary updates to capture the rotor wake effects. Once this baseline work 
is complete, full hover rotor performance map correlations should be completed for a variety of rotors 
in order to substantiate the promise of computational improvements due to the ZAPR3D coupling con- 
cept. 
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4.0 CONCLUSIONS 


The primary objective of this work was to demonstrate the feasibility of a new potential/viscous 
flow coupling procedure for reducing computational effort while maintaining solution accuracy. The 
success of such a procedure should allow for more exact viscous simulations for simple geometries on 
a workstation environment and full helicopter configurations on current super computers. In this zonal 
procedure, the flow field is partitioned into zones which utilize approximate equations appropriate to 
the flow physics within each zone. The coupling concept is based on the premise that any computa- 
tional method can only produce valid results within the approximations of the physical model employed 
in its construction. Therefore, the interfacing boundary surface between the potential and viscous flow 
code domains must be very nearly potential in order to achieve a successful (accurate) result. This is 
accomplished in the current technique by requiring that beyond the first iteration the potential flow 
surface panel strengths are obtained from the NS solution at a smooth inner fluid boundary. These 
fluid surface panel values are then used to compute the outer velocity boundary conditions for the NS 
calculation. The iteration between the potential flow solution and the NS solution continues in a closed 
loop until flow convergence, or allowable iteration limits, are obtained. 

This closed-loop, overlapped, velocity-coupling concept has been developed in a new two-dimen- 
sional code, ZAP2D (Zonal Aerodynamics Program-2D), a three-dimensional code for wing analysis, 
ZAP3D (Zonal Aerodynamics Program-3D), and a three-dimensional code for isolated helicopter rotors 
in hover, ZAPR3D (Zonal Aerodynamics Program for Rotors-3D). ZAP2D couples a simple airfoil 
panel code, POT2D, with the airfoil Navier-Stokes code, ARC2D. Similarly, ZAP3D couples a modi- 
fied version of VSAERO with ARC3D. Finally, ZAPR3D, couples the modified version of VSAERO 
with the lifting surface HOVER code and the ARC3D code that has been modified to include rotational 
effects. The current status and findings for each of these programs are briefly summarized below. 

Detailed results developed and obtained under this contract for ZAP2D applied to a NACA 0012 
airfoil at low angle of attack was reported in Reference 24. In this report, studies for a range of angles 
of attack and Mach number are thoroughly discussed. Computational effects of grid domain sizes from 
25 to 0.10 chord lengths have been investigated for a NACA 0012 airfoil with ZAP2D. Comparisons 
with large domain ARC2D solutions and with experimental data have shown that the required domain 
size can be reduced to a few tenths of a percent chord for the low Mach and low angle of attack cases 
and to less than 2-5 chords for the high Mach and high angle of attack cases while maintaining solution 
accuracies to within a few percent. The smaller domain sizes result in a reduction in the required 
number of grid points and also a reduction of the number of required iterations for a converged solu- 
tion; consequently, at this stage of the development, ZAP2D has demonstrated CPU reductions by 
factors of 2-4 compared with ARC2D. 

Future enhancements to ZAP2D should include an automated procedure to determine the opti- 
mum location and shape of the inner and outer coupling boundaries. Such improvements in ZAP2D 
should routinely provide a reduction in CPU requirements for all cases by a factor on the order of 4 
compared with the N-alone simulation. 

Detailed calculations by ZAP3D and ARC3D have been investigated for a rectangular planform 
wing of aspect ratio 5 with NACA 0012 airfoil sections. All calculations included a ffeestream angle 
of attack of 8.75°, a Mach number of 0.12, and a Reynold’s number of 1.5 million. In this work, grid 
domain sizes from 17 chord lengths to 0.41 chord lengths have been studied. Converged values for 
C L , C D , and C^ were obtained from the large domain (R 0 =17) ARC3D calculation and served as 
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target values for the ZAP3D simulations for reduced grid domain sizes. Current calculations show that 
the grid domain size for ZAP3D can be reduced to 0.7 chord lengths with less than a 3% error in C L 
compared with the large domain ARC3D result. Additional reductions in the required computational 
domain for ZAP3D are expected as the method is further developed and refined. Nevertheless, the 
current ZAP3D calculation for an outer domain radius of about 1.2 chords represents a speed-up in 
CPU time over the ARC3D large domain calculation by about a factor of 2.5. This improvement is 
achieved for less than a 0.5% deviation in C L , 10 counts change in C D , and 0.0015 variation in C^. 

Future work with ZAP3D should address modeling improvements in the downstream wake region 
in order to further reduce the domain size requirement and improve the computational efficiency of the 
code. Higher angle of attack and Mach number cases also need to be investigated before the overall 
effectiveness of the zonal coupling procedure utilized in ZAP3D can be validated. 

Preliminary hover calculations with ZAPR3D have been completed for a two-bladed, rotor with a 
NACA0012 airfoil and an aspect ratio of 5. As noted above, ZAPR3D combines methodology from 
ZAP3D (potential flow updates at the inner boundary surface, S ■), HOVER (lifting surface hover meth- 
od to initialize the rotor solution and to provide the wake and secondary blade influences that are out- 
side of S;), and a modified ARC3D (viscous solution within the outer domain, S„). The addition of the 
rotational flow terms to ARC3D were validated by comparison with experimental data for the non- 
lifting hover case. As shown in this report, computed and experimental chordwise pressure distribu- 
tions are comparable for the non-lifting case. Coupled ZAPR3D calculations for a lifting case with a 
collective setting of 8 degrees have included selected domain sizes from approximately 10 chord down 
to 1 chord lengths. While a full 180 degree azimuthal grid calculation with periodic boundary condi- 
tions was not able to capture and maintain the full wake effects and did not achieve a favorable com- 
parison with experimental pressure data, a ZAPR3D simulation for a fine grid restricted to a reduced 
domain of about 2 chord lengths was able to capture the wake effects and did compare accurately with 
the experimental pressure data. Current preliminary comparisons of the radial thrust loading computed 
by HOVER and various ZAPR3D calculations show good agreement in the maximum thrust loading 
due to the tip vortex passage but indicate that this peak value is located further outboard in ZAPR3D 
than in HOVER. Furthermore, the ZAPR3D loading near the root also indicates that the boundary 
conditions imposed there may not be consistent with the HOVER simulation. Comparisons of the 
residual history for a ZAPR3D simulation with no inside loop updates (ITI=0) and a ZAPR3D simula- 
tion with inside updates every 250 iterations (ITI = 7) illustrate that the inner surface updates are cor- 
rectly modeled but also indicates that the domain size needs to be reduced considerably in order to 
demonstrate the convergence improvement that is found in the ZAP3D wing simulations. Until this 
additional work is complete, estimations of any reductions in computational effort cannot be made. 

Future ZAPR3D work should concentrate on the detailed validation of the effectiveness of the 
inner boundary updates to capture the rotor wake effects and to systematically reduce the required 
outer boundary domain size. Once this baseline work is complete, full hover performance map corre- 
lations should be completed for a variety of rotors in order to substantiate the promise of computational 
improvements due to the ZAPR3D coupling concept. 
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Fig. 1. Zonal Representation of the Physical Flow Field. 



Fig. 2. Zonal Methodology Flow Chart. 
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Fig. 5. Rotor Blade Vortex-Lattice Model. 


Fig. 6. Global Hover Wake Model. 
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Fig. 9. Computed Force and Moment Variation vs. Angle of Attack for NACA 0012 Airfoil with 
Domain Size Variation; M=0.3, R e =3.0 Million. 






(a) Lift vs. Drag Coefficient 



Fig. 10. Drag Polar and Pitching Moment Curve with Domain Size Variation; M=0.3 
and R e =3.0 Million. 
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Fig. 11. 
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(b) ZAP2D with Grid F 


Fig. 12. Comparison of ZAP2D and ARC2D Computed Mach Contours;M=0.3, R e =3.0 Million 
and a = 12.86°. 
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Comparison of ZAP2D and ARC2D Computed Pressure Contours; M=0.3, R e =3.0 Million 
and a= 12.86°. 
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Fig. 14. Streamlines around the NACA 0012 Airfoil; M=0.3, R e — 3.0 Million and a— 12.86°. 





Fig. 15. Comparison of ZAP2D (R^O.17) and ARC2D (R„=25) Computed Loads and Residual 
Iteration History; M=0.3, R o =3.0 Million and a= 12.86°. 
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(c) ZAP2D (Grids B-H) and 
ARC2D (Grid A) Computed 
Moment Coefficient 


Fig. 16. Computed Force and Moment Variation with Domain Size for a NACA 0012 Airfoil at 
Various Mach Numbers and at Angles of Attack near C,^^. 
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Fig. 17. Error i; 
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Fig. 18. Comparison of ZAP2D (Grids C and G), ARC2D (Grid A) and Experimental Force and 
Moment Variation with Mach Number. 






Fig. 19. Comparison of ZAP2D and ARC2D Computed Mach Contours for a NACA 0012 Airfoil at 
Mach 0.7, a=4.75°. 
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Fig. 21. Comparison of ZAP2D (Grid C) and ARC2D (Grid A) Computed C f Iteration History for a 
NACA 0012 Airfoil at Mach 0.7, a=4.75°. 



Fig. 22. Comparison of ZAP2D (Grid C) and ARC2D (Grid A) Solution Residual History for a 
NACA 0012 Airfoil at Mach 0.7, a=4.75°. 
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(a) Horizontal Component 



(b) Vertical Component 


Fig. 23. Baseline ZAP3D and ARC3D Computed Velocity Component Normalized by Local Speed 
of Sound for an AR 5 Rectangular Wing at Mach 0.12, a=8.75 and R e 1.5 Million. 









Fig. 24. ZAP3D and ARC3D Computed Velocity Component Normalized by Local Speed of Sound 
with Short S; for an AR 5 Rectangular Wing at Mach 0.12, a=8.75 and R e 1.5 Million. 
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(c) Moment Coefficient 



(b) Drag Coefficient 



(a) Lift Coefficient 


Fig. 25. ZAP3D and ARC3D Computed Force and Moment Variation with Domain Size for an AR 
5 Rectangular Wing at Mach 0.12, a=8.75 and R e 1.5 Million. 
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Fig. 27. Computed ai 
1.5 Million. 
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Fig. 28. History 
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Fig. 29. 
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Fig. 31. Computational Model of the Wing/Vortex Validation Case. 



(a) Lift Induced by a Vortex -generating Wing 



(b) Lift Induced by a Prescribed Vortex Strength NACA 0012 Wing for a Non-lifting Rotor Case 


Fig. 32. Comparisons of Calculated lift on a NACA 0012 Wing. 
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Fig. 33. A Partial View of the Grid Generated around a NACA 0012 Wing for a Non-lifting Rotor Case. 
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Fig. 34. Calculated Pressure Coefficients Based Upon a Local Dynamic Pressure along the Span Direction (M=0.52, R e =2.32 Million, 
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(b) The Span Location Y=0.68 
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(c) The Span Location Y=0.80 



(d) The Span Location Y=0.89 
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(e) The Span Location Y=0.96 



(a) Axial Coordinates 



Fig. 35. Comparison of Experimental and Calculated Tip-vortex Geometry for the Un-twisted Rotor 
of Ref. 10. 
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Fig. 36. Final Calculated Relaxed Wake Geometry at Constant Azimuth for the Un-twisted Rotor of 
Ref. 10. 



Fig. 37. Comparison of Calculated and Experimental Hover Performance for the Un-twisted Rotor of 
Ref. 10. 
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(b) Flow Field Outside of Sj 


Fig. 38. The Downwash Predicted by HOVER for the Viscous Calculations at an Upstream Plane. 
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Fig. 39. Comparison of the Calculated and Experimental Chorwise Pressure Distribution at r/R=0.68 for a Rotor in Hover (M=0.44, 
Re= 1.92 Million, 0=8°). 




Fig. 40. Selected Planes and Grid for a Reduced Domain (Zonl) Simulation. 



Fig. 41. Comparison of the Calculated Radial Thrust Distribution by HOVER and ZAPR3D 
(Zonl) (M=0.44, Re=1.92 Million, 0 = 8 °). 
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Fig. 44. 


Convergenc e History for a Two-Bladed Rotor (M=0.44, Re=1.92 



Fig. 45. Comparison of HOVER and ZAPR3D Computed Radial Thrust Distribution for a 
Two-Bladed Rotor (M=0.44, Re=1.92 Million, 0=8°). 
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Fig. 46. Comparison of ZAPR3D Calculated and Experimental Chordwise Pressure Distribution for a Two-Bladed Rotor (M=0.44, 

Re= 1.92 Million, 0=8°). 
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Fig. 46. Continued. 
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Fig. 46. Continued. 
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Fig. 46. Continued. 
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Fig. 46. Concluded. 
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